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Abstract Geographic variation of morphology is an 
important topic of evolutionary biology, and research 
on geographic variation can provide insights on the 
formation, evolution, and adaptation of species and 
subspecies. The vertebrate skull is a developmentally 
and functionally complex morphological structure with 
multiple functions, that is susceptible to vary according 
to selection pressure. In this study, geographic variations 
in skull morphology of Batrachuperus karlschmidti from 
four different geographic populations (Shade, Gexi, 
Shangluokema, and Xinduqiao) were examined via 
geometric morphometrics. No significant differences 
were found among these populations with regard to 
skull size; however, significant variation was found 
in skull shape. The most notable shape changes are 
the relative sizes and positions of the frontal, maxilla, 
pterygoid, and vomer. Skull shape changes were not 
related to allometry. However, due to limitation of 
sample populations and size, the results of this study 
need to be further verified by more sample populations 
and individuals in the future. The results of this study 
contribute to our knowledge about these aspects of 
morphological variability in this species as well as in 
hynobiid salamanders. 
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1. Introduction 


Geographic variation and the variation among populations in 
genetically based traits across the natural geographic range of 
a species (Mary, 1977; Olvido and Mousseau, 2012) has been a 
central theme in evolutionary biology (Fonseca and Asttia, 2015). 
The traits of geographic variation in amphibians include many 
aspects and have been widely studied, such as morphology 
(eg, Valenzuela-Sanchez et al., 2015; Oromi et al., 2016; Yang 
et al, 2017; Liu et al, 2018; Amado et al., 2019), behavior (eg., 
Rodriguez-Tejeda et al, 2014; Wise and Jaeger, 2016), physiology 
(e.g, Riddell and Sears, 2015; Twomey et al., 2015; Barria and 
Bacigalupe, 2017; Xiong et al, 2018), and life history (e.g, Kopp 
and Baur, 2000; Miaud et al, 2001; Morrison and Hero, 2003; 
Oromi et al., 2014, 2016; Liao et al., 2016; Wang et al., 2019). 
Geographic variation in morphology is the most widespread 
phenomenon (Avise, 2000), which reflects phenotypic responses 
to environmental gradients, genetic variation and evolutionary 
history of populations and species and may indicate local or 
regional changes in environmental conditions (Ellison et al., 
2004). Genetic components and environmental conditions 
were considered to be the main reasons resulting in geographic 
variation (Olvido and Mousseau, 2012). Studying geographic 
variation can provide insights on the formation and evolution 
of species and subspecies, the generation of biodiversity, the 
adaptive mechanisms of species, as well as the predication of the 
adaption of species to future environmental changes (Losos and 
Glor, 2003; Ellison et al., 2004; Olvido and Mousseau, 2012). 

The vertebrate skull, a developmentally and functionally 
complex morphological structure with multiple functions 
such as perception, protection, competition, reproduction, and 
predation (Hanken and Hall, 1993; Ivanović et al., 2013), is a 
target of selection (Alarcén-Rios et al., 2017). The morphology 
of the vertebrate skull easily allows for variations among 
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populations of a species due to differences of environmental 
conditions. Geometric morphometrics (Rohlf and Marcus, 
1993) has been widely and successfully applied to research 
morphological variation of the vertebrate skull, which is 
particularly true for Urodela. Ivanović et al. (2009) showed 
significant variation in skull size among populations of 
Mesotriton alpestris. Ivanović et al. (2012) found statistically 
significant variation in skull size and skull shape among 
populations of Lissotriton vulgaris. Ivanović et al. (2013) reported 
significant divergence in skull shape but not in skull size 
between mitochondrial DNA lineages of Triturus karelinii. Üzüm 
et al. (2015) found statistically significant variation in skull size 
and shape among populations of Ommatotriton ophryticus and O. 
vittatus, and the variation in skull shape in these two species was 
almost entirely due to size-dependent, allometric shape changes. 
However, no research of geometric morphometrics on skull of 
hynobiid salamanders has been reported to date. 

The Chiala mountain salamander, Batrachuperus karlschmidti, 
was named by Liu (1950) based on specimens from Chiala 
[= Jiana], Luhohsien [= Luohuo County], Sikang [= Sichuan 
Province], China. It is differentiated from all other congeneric 
species (B. pinchonii, B. londongensis, B. tibetanus, and B. 
yenyuanensis) mainly by the uniformity of its dark olive-colored 
body, with a total lack of spots and marblings (Liu, 1950). 
Several scholars have considered B. karlschmidti to be a junior 
synonym of B. tibetanus (Fei et al., 1983; Fei et al, 2006); however, 
the morphology, development, and habitat preferences (Zhao 
and Jiang, 1988), cytogenetics (Yang, 1992), and molecular 
biology (Fu et al., 2001; Li et al, 2004; Fu and Zeng, 2008; Lu et al., 
2012; Xia et al, 2012) support its validity. This species is widely 
distributed in the mountain streams of Sichuan Province (Fu 
and Zeng, 2008; Lu et al., 2012; Xia et al, 2012) and Qinghai 
Province (Ding et al, 2014). This species is adapted to cold water 
and can be found under stones of small mountain streams from 
1500 m to 4250 m, and the breeding season lasts from May 
to the first part of August (Liu, 1950). This species is a poorly 
known and understudied species, and have been studied rarely. 
In this study, geographic variation in the skull morphology 
among four populations (Shade, Gexi, Shangluokema, and 
Xinduqiao populations) of B. karlschmidti was analyzed using 
geometric morphometric methodology. The following questions 
were preliminary tested: (i) whether geographic divergence 
existed in skull size and shape, and (ii) whether divergences in 
skull shape are related to the skull size. This study improves our 
understanding of the evolution of skull shape and processes 
that lead to morphological diversification. 


2. Materials and methods 


2.1. Samples and skull preparation A total of 24 skulls of 
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B. karlschmidti were used in this study. Species were collected 
from four populations of Sichuan Province, China (Figure. 1): 
Shade (29.5123°N, 101.4374°E, altitude 3447 m), Gexi (30.9696°N, 
101.2229°E, altitude 3622 m), Shangluokema (31.6168°N, 
100.7780°E, altitude 3644 m), and Xinduqiao (29.7906°N, 
101.5623°E, altitude 3743 m). Each population included six 
specimens (three males and three females). Only adult specimens 
(identified via examination of the gonad development) were 
included in this analysis to avoid additional ontogenetic 
variation. Upon arrival at the laboratory, animals were 
euthanized via submergence in a buffered MS-222 solution 
and were then stored in 10% formalin. Specimens were sexed by 
inspection of the gonads. The skulls were prepared by clearing 
and double-staining using a bone and cartilage staining 
procedure (Hanken and Wassersug, 1981), which has been 
briefly summarized by Xiong et al. (2013) and Xiong et al. (2016). 
The cleared and stained skulls were stored in glycerol and 
deposited in the Henan University of Science and Technology 
Museum (HNUSTM). 

All applicable international, national, and/or institutional 
guidelines for the care and use of animals were strictly 
followed. All animal sample collection protocols complied with 
the current laws of China. 


2.2. Data acquisition Pictures of the stained skull specimens 
were taken via a mounted OLYMPUS DP26 digital camera. 
The skulls were positioned in the center of the optical field to 
reduce and equalize distortion (Cvijanovié et al, 2017). TpsUtil 
ver. 1.74 (Rohlf, 2017) was employed to create TPS files that 
contained all coordinate datasets. Landmarks, which enabled 
the definition of shapes, were digitized on the pictures of each 
skull specimen using TpsDig2 ver. 2.16 (Rohlf, 2010). Nineteen 
landmarks were used for the dorsal and 21 landmarks for the 
ventral view of the skull (Figure. 2). All landmarks were placed 
on the left side of the skull, and were digitized three times by 
the same person (QQL). Digitizing hemi-skull was not only 
done to avoid redundancy information in symmetric structures, 
but also to reduce the number of analyzed parameters in the 
statistical processing of the data set (Abramon et al., 2017). In 
addition, although directional asymmetry and fluctuating 
asymmetry of shape are widespread in the animal kingdom, 
these are very small (Klingenberg, 2015). Thus, the potential 
left-right differences, directional asymmetry and fluctuating 
asymmetry were ignored in this investigation. The definition of 
each landmark is listed in Table 1. 


2.3. Geometric morphometrics A full Procrustes fit in 
MorphoJ v 1.06d was performed to eliminate effects of 
orientation, position, and size (Klingenberg, 2011), which 
provided two new set of variables: Centroid size (CS, the square 
root of the sum of squared distance of all landmarks from 
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Figure 1 Localities of samples used in the present study. 


Figure 2 Landmarks digitized on the skull in dorsal (A) and ventral (B) views. 
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their centroid) values (Dryden and Mardia, 1998) and a matrix 
of shape coordinates (Procrustes coordinates, representing the 
difference between the consensus specimen and each sample). 
The centroid size value was computed for each specimen. 
Procrustes ANOVA (Klingenberg and McIntyre, 1998; 
Klingenberg et al, 2002) was used to estimate the measurement 
error. When the measurement error was not significant, 
the mean values and co-ordinate of each individual were 
calculated and used for the following analyses. The differences 
of CS among populations were tested by analysis of variance 
(ANOVA) using the RRPP package (Version 0.5.1) in R (Version 
3.6.2). Principal component analysis (PCA) was performed on 
the covariance matrix of the Procrustes coordinates to explore 
and visualize the main trends of variation among populations. 
To further quantify the shape differences between populations, 
Procrustes distances, which is an absolute measure of the 
magnitude of shape difference between two configurations, 
were calculated for each pair of populations. The statistical 
significance of differences in the mean shape between 
populations was estimated by permutation test based on 10 000 
iterations. Multivariate regression of Procrustes coordinates 
on log CS was conducted to estimate the relationship between 
shape and size (allometry) with 10 000 permutation among and 
within populations (Klingenberg, 2011). 

Due to the limited specimen number, differences of size and 
shape between sexes were not explored in this analysis. The 
analyses were run separately for the dorsal and ventral data. 
All geometric morphometric analyses and visualization were 
conducted with MorphoJ v 1.06d (Klingenberg, 2011), unless 


stated otherwise. 
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3. Results 


3.1. Measurement error Measurement error explained the 
smaller portion on both dorsal and ventral size and shape 
variation (dorsal size 0.1%, dorsal shape 2.0%, ventral size 0.1%, 
and ventral shape 2.9%) (Table 2). Thus, the digitalizing error 
was not considered in our data, and the mean values of three 
times of each individual were calculated and used for the 


following analyses. 


3.2. Skull size variation The mean CS values of Shade and 
Xinduqiao populations were similar, and larger than those of 
Gexi and Shangluokema populations, which were similar, in 
both dorsal and ventral skull (Table 3). Although the mean CS 
values showed variations, the analysis results of ANOVA with 
permutation showed no significant variation neither for dorsal 
(F = 03707, P = 0.539) nor ventral skulls (F = 0.6526, P = 0.419) 


among populations. 


3.3. Skull shape variation The results of Procrustes ANOVA 
showed that the shape differed significantly among populations 
for both dorsal (F = 3.58, P < 0.0001) and ventral (F = 2.76, 
P < 0.0001) skulls. For the dorsal skull, the first three principal 
components described 72.56% of the total shape differences. 
PCI (41.64%) showed the shape change of frontal (landmark 2) 
(Figure 3); PC2 (23.38%) mainly reflected the shape changes of 
squamosal (landmarks 8, 9, 10) and exoccipital (landmarks 4, 
5, 6) (Figure 3); and PC3 (7.54%) represented the shape changes 
of squamosal, exoccipital and processus dorsalis praemaxillaris 
(landmark 11) (Figure 4). In the plot of the first two axes, 
Shade population largely overlapped with Shangluokema 
and Xindugiao populations, but could be clearly distinguished 


Table 1 Definitions of the landmarks digitized on the dorsal and ventral skull and their anatomical descriptions. 


Dorsal skull Ventral skull 
Number Description Number Description 
1 Tip of the snout 1 Tip of snout 
2 Posterior point of frontal 2 Anterior point of suture between vomer and premaxilla 
3 Posterior point of parietal 3 Posterior point of vomer 
4 Rightmost point of exoccipital 4 Posterior point of parasphenoid 
5 Posterior point of exoccipital 5 Posterior point of operculum 
6 Leftmost point of exoccipital 6 Anterior point of operculum 
7 Outside point of parietal 7 Outside point of operculum 
8 Posterior point of squamosal 8 Outward point of pterygoid 
9 Anterior point of squamosal 9 Anterior point of pterygoid 
10 Outside point of squamosal 10 Inside point of pterygoid 
11 Posterior point of the processus dorsalis praemaxillaris 11 Posterior point of suture between vomer and maxilla 
12 Posterior point of maxilla 


13 Suture between maxilla and premaxilla 
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from Gexi population; Gexi population overlapped with 
Shangluokema population, but could be clearly distinguished 
from Shade and Xindugiao populations. Shangluokema 
population only could be clearly distinguished from Xinduqiao 
population (Figure 3). The ordination of PC1 vs. PC3 enabled 
clear distinction of Shade and Xindugiao populations from 
Gexi and Shangluokema populations, respectively; however, 
Gexi and Shangluokema populations largely overlapped 
(Figure 4). The pairwise Procrustes distances (Table 4) indicated 
that the distance between Gexi and Shangluokema populations 
was smallest, while the distance between Gexi and Xinduqiao 
largest. Compared with the overall mean configuration, Shade 
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population had slight backward-extending frontal (Figure 5A); 
Gexi population had backward-extending frontal (Figure 5B); 
Shangluokema population had slight laterally extending frontal 
(Figure 5C); Xindugiao population had forward-extending 
frontal (Figure 5D). 

With regard to the ventral skull, the first three principal 
components described 59.047% of the total shape differences. 
PCI (30.60%) reflected the shape changes of maxilla, pterygoid 
and vomer (Figure 6); PC2 (17.45%) represented shape changes 
of pterygoid, maxilla, and vomer (Figure 6); PC3 (11.00%) mainly 
reflected shape changes of petergoid, vomer, and maxilla 
(Figure 7). In the plot of the first two axes, except for Shade and 


Table 2 Procrustes ANOVA of skull size and shape in Batrachuperus karlschmidti from three replicates. 


Effect Sums of squares Mean squares df F P 
Dorsal size 

Individual 57.781526 2.51224 23 1069.06 <0.0001 
Residual (digitization) 0.112797 0.00235 48 

Dorsal shape 

Individual 0.10414362 0.0002515546 414 49.77 <0.0001 
Residual (digitization) 0.00436735 0.0000050548 864 

Ventral size 

Individual 83.717186 3.639878 23 213.38 <0.0001 
Residual (digitization) 0.081704 0.001702 48 

Ventral shape 

Individual 0.1071064 0.0002116727 506 34.7 <0.0001 
Residual (digitization) 0.00644202 0.000006 1004 1056 
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Figure 3 Shape variations in the dorsal skull view. A. Scatter plot of PC1 versus PC2; B-E. Wireframe deformations depict the shape 
changes along the respective PC axes. 
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Xinduqiao populations, which largely overlapped and could not 
be distinguished, Gexi and Shangluokema populations could be 
well distinguished and could also be distinguished from Shade 
and Xinduqiao populations (Figure 6). The ordination of PC1 
vs. PC3 showed that Shangluokema and Gexi populations 
largely overlapped, and Shade and Xindugiao populations also 
largely overlapped; however, Shangluokema population could 
be clearly distinguished from Shade and Xindugqiao populations 
(Figure 7). The pairwise Procrustes distances indicated that 
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the distance between Shade and Xinduqiao populations was 
smallest, and the distance between Shade and Shangluokema 
population was largest (Table 4). Compared with the overall 
mean configuration, Shade population had backward- 
extending vomer (landmarks 2, 3, 11), laterally extending 
maxilla (landmarks 12, 13) and pterygoid (landmarks 8, 9, 10) 
(Figure 8A); Gexi population had longer maxilla (Figure 8B); 
Shangluokema population had shorter maxilla and adducent 
pterygoid (Figure 8C); Xindugiao population had adducent 


Table 3 Descriptive statistics of dorsal and ventral skull sizes (means + standard error) of Batrachuperus karlschmidti among different 
populations. Skull sizes are given as centroid sizes (CS); N: number of specimens. 


Dorsal Ventral 
Populations N 
Mean+SE Ranges Mean+SE Ranges 
Shade 6 16.498+0.429 14.903-18.018 18.588+0.500 16.680-20.425 
Gexi 6 15.503+0.468 14.686-17.609 17.295+0.533 16.220-19.733 
Shangluokema 6 15.427+0.160 14.996-15.946 17.289+0.285 16.368-18.383 
Xindugiao 6 16.180+40.223 15.342-16.654 18.044+0.273 17.104-18.609 


Table 4 Procrustes distance of dorsal and ventral skull between pairwise geographical populations and its respective P-values after 
permutation test with 10 000 permutation rounds. Significant of Procrustes distances are shown in bold. 


Procrustes dorsal distances 


Procrustes ventral distances 


populations 


Gexi Shangluokema Shade Gexi Shangluokema Shade 
Shangluokema 0.0196 0.0301 
Shade 0.0338 0.0283 0.0378 0.0453 
Xinduqiao 0.0551 0.0422 0.0306 0.0296 0.0367 0.0194 
D oA 
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Figure 4 Shape variations in the dorsal skull view. A. Scatter plot of PC1 versus PC3; B-E. Wireframe deformations depict the shape 


changes along the respective PC axes. 
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Figure 5 Wireframe deformations depict the shape changes of the overall mean versus population means in dorsal skull. A. Shade 
population mean versus overall mean; B. Gexi population mean versus overall mean; C. Shangluokema population mean versus 
overall mean; D. Xindugiao population mean versus overall mean. 
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Figure 6 Shape variations in the ventral skull view. A. Scatter plot of PC1 versus PC2; B-E. Wireframe deformations depict the shape 
changes along the respective PC axes. 


pterygoid (Figure 8D). ventral: total sums of squares = 0.0252, Predicted = 6.3483%), 
but shape changes did not significantly correlate with changes 
in size both in dorsal (P = 0.2103) and ventral skull (P = 0.1190). 
Within each population, multivariate regression showed that 
showed that size high percentages of prediction were obtained shape changes significantly correlated with changes in size 


(dorsal: total sums of squares = 0.0226, Predicted = 5.8905%; only in dorsal skull of Gexi population (P = 0.0193), but not 


3.4. Static allometry Multivariate regression of the Procruste 
coordinates on log CS (pooled regression within populations) 
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significantly in dorsal and ventral skull of other populations 
(P > 0.05). 


4. Discussion 


To the best of our knowledge, this study is the first to use 
geometric morphometric methods to describe size and shape 


0.024 


0.007 


PC3 (11.00%) 
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variations of the skulls of hynobiid salamanders. Although the 
sample populations and size are limited, the presented results 
shed new insight into important issues of geographic variation 
in hynobiid salamanders. The outcome of the present study 
demonstrates that geographic variation occurs in skull shape 
but not in skull size. The skull shape changes did not result from 


allometry. 
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Figure 7 Shape variations in the ventral skull view. A. Scatter plot of PC1 versus PC3; B-E. Wireframe deformations depict the shape 


changes along the respective PC axes. 


Figure 8 Wireframe deformations depict the shape changes of the overall mean versus population means in ventral skull. A. Shade 
population mean versus overall mean; B. Gexi population mean versus overall mean; C. Shangluokema population mean versus 
overall mean; D. Xindugiao population mean versus overall mean. 
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Centroid size was used as an indicator of the overall 
body-size of specimens in geometric morphometric analyses 
(Zelditch et al., 2012). Geographic variation of skull size 
(centroid size) varies from species to species. For example, 
geographic variations of skull size are present in Triturus 
cristatus superspecies (Ivanović et al., 2008), Mesotriton alpestris 
(Ivanović et al, 2009), Lissotriton vulgaris (Ivanović et al., 2012), 
Ommatotriton ophr yticus and O. vittatus (Üzüm et al, 2015), but are 
absent in Triturus karelinii lineages (Ivanović et al, 2013). In this 
study, the centroid size of both dorsal and ventral skull did not 
significantly vary among four B. karlschmidti populations. In 
taxa with indeterminate growth, age, environmental conditions, 
and habitat resources largely affect body size (Arntzen, 2000; 
Gotthard, 2001; Cogalniceanu and Miaud, 2003; Abramov et al., 
2017). In this study, only adult individuals were analyzed, and 
habitat resources of four populations were similarly based on 
our field survey. Although the environmental conditions (e.g, 
altitude, longitude, and latitude) were different, the differences 
were small. These result indicate a lack of variation of skull size 
among the four tested B. karlschmidti populations. 

Since the vertebrate head is a structure with relevance 
in ecological and social functions, the shape of the skull is 
susceptible to selection pressures (Alarcén-Rios et al., 2017). 
Previous analyses on variations of the urodela skull shape 
showed a high level of variation in skull shape between 
populations, and notable shape changes are varied. For 
example, the most prominent shape changes were in the skull 
base, palatal, snout, and vomeral teeth in Triturus karelinii 
(Ivanović et al, 2013), in the skull base, palatal, vomer, snout, 
maxilla, and praemaxillae in Lissotriton vulgaris (Ivanović et al., 
2012), as well as the exoccipital, snout, squamosal, and frontal 
in Mesotriton alpestris (Ivanović et al, 2009). In this study, the 
most notable shape changes were found in the frontal, maxilla, 
pterygoid, and vomer among the four tested populations. 
These shape changes can clearly distinguish most populations. 
Compared with dorsal skull, ventral skull showed higher 
variation, this is consistent with the findings of other studies 
(eg, Ivanović et al, 2011, 2013; Üzüm et al, 2015). The variation 
of skull shape within species may result from static allometry, 
genetic differentiation, and ecological preferences (Ivanović 
et al, 2009, 2012, 2013; Ivanović and Arntzen, 2014; Üzüm et 
al., 2015). The results of multivariate regression showed no 
significant divergence in allometry of dorsal and ventral 
skull among populations, and regression within population 
showed that significant divergence in allometry was only 
found in the dorsal skull of the Gexi population. These results 
suggest that allometry is not the influencing factor in the skull 
morphological variation of B. karlschmidti. The genetic data and 
ecological preferences (e.g, dietary) of the individuals of these 
four B. karlschmidti populations are absent. Thus, the variation 
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of skull shape among populations of B. karlschmidti cannot be 
explained by existing data, and more information of these 
populations is needed. 

Research on geographic variation is usually based on a 
larger number of populations and individuals (e.g, Ivanović et 
al., 2008, 2009, 2013; Üzüm et al., 2015). Here, the skull size and 
shape variation of B. karlschmidti was studied based on limited 
sample populations and size. The obtained results of this study 
need to be verified/explored by more sample populations and 
larger sample sizes in the future. Geographic variation may 
result from genetic components and environmental conditions 
(Olvido and Mousseau, 2012). In this study, the genetic data 
of four populations are absent. Such data can not only be 
used to analyze the mechanism of variation but also provides 
phylogeographic structuring to analyze the phylogenetic signal 
of variation. Then, genetic studies should also be included in the 
following research. 


5. Conclusion 


In conclusion, geographic variation does not affect skull size, 
which may result from a lack of differences of affecting 
factors (age, environmental conditions and habitat resources), 
but rather affects skull shape. Due to the limitation of sample 
populations and their size in this study, to better understand 
the mechanism and skull shape disparity pattern, considerable 
efforts should be directed toward a broader range of 
populations and individuals, and the genetic diversity of this 
species should be studied in the future. 
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